Method for determining physical properties of a multilayered periodic structure

ABSTRACT

There is provided a method of calculating physical properties of a periodic structure. At least one physical property related to reflectivity or transmittance of a periodic structure is measured, and then, at least one physical property related to reflectivity or transmittance of a virtual periodic structure is calculated to obtain corresponding physical properties from the virtual periodic structure. The at least one calculated physical property is compared with the at least one measured physical property. When the virtual periodic structure is horizontally divided into a plurality of layers, at least three substances can have horizontally repeated periods in the middle layers of the divided structure. In accordance with an embodiment of the present invention, the microscopic formation of the periodic structure including a native oxide layer formed on the periodic structure or an intentionally formed surface coating layer thereon can be nondestructively and accurately tested.

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of Korean Patent Application Nos.10-2007-0052806 filed on May 30, 2007; 10-2007-0077351 filed on Aug. 1,2007; 10-2007-0104495 filed on Oct. 17, 2007; and 10-2007-0104496 filedon Oct. 17, 2007, the disclosures of which are hereby incorporatedherein by reference.

BACKGROUND OF THE INVENTION

1. Technical Field

Embodiments of the present invention described herein relate to a methodof testing a periodic structure, e.g., to determine a shape of astructure having a feature repeated in one or more dimensions, themethod performed such as through non-destructive testing via measurementof reflectivity or transmittance.

2. Discussion of Related Art

Generally, to fabricate electronic devices, such as semiconductordevices or display devices, processes of cleaning, thin-film growing,photolithography and thin-film etching are repeated many times toproduce the consumer products. For example, in the photolithographyprocess, a circuit of a mask where an image to be fabricated is formedand is transferred to a photosensitive material (photo resist) to form apattern, and the pattern is used as an etch barrier to form a desiredcircuit on a thin film.

In the semiconductor and display devices fabricated by using thephotolithography process, the desired circuit needs to be transferred tothe thin film in an accurate shape in each step. This is possible basedon the accuracy of the photolithography process. That is, only when theshape of a desired pattern is accurately transferred to a photo resist,and the resist layer properly functions as the etch barrier, an accuratecircuit can be formed on the thin film. That is, the accurate pattern isto be formed by the photo resist before the circuit is formed on thethin film, and this can be confirmed by a testing process.

To test a pattern, there has been generally used a method of opticallyobserving a shape of a semiconductor device using a pattern tester, forexample. However, since the resolution of the pattern tester can beinsufficient for determining the shapes of “nano-level” patterns whichmeasure only a few nanometers in length. Using a pattern tester, it isdifficult to perform an accurate analysis. To solve such a drawback, inthe semiconductor research and production line, there has been used amethod of analyzing a specific shape using equipment such as an electronmicroscope.

However, when an electron microscope is used, since a section of asemiconductor device is to be cut for the analysis of a shape thereof,the semiconductor device as fabricated cannot be used again. Moreover,since the measurement is to be conducted under vacuum environment, itcan take an excessively long time to obtain a result of the measurement.It may also be impossible to select particular regions of a sample to bemeasured. Due to the aforementioned drawbacks, the electron microscopehas a limit in its practical use in the production line.

To address the aforementioned drawbacks, the technology using an opticalmeasurement method has been developed and includes, for example, amethod of using an approximate expression called the Effective MediumApproximation (EMA). A calculation method using the EMA has the problemin that, since an approximation is obtained by only a volume ratio ofconstituent substances in a given period, regardless of a detailed shapeof a structure, it never distinguishes the detailed shape of thestructure. That is, since the shape of each pattern of a circuit with aperiodic structure is not specifically distinguished and only the volumeratio of constituent substances in a given period is distinguished, thedifference between the real structure and the measured structure issignificant. Specifically, in the periodic structure, since thecalculation method using the EMA cannot clarify the different periodicstructures if their volume ratios are same, a new optical measurementmethod is really needed.

SUMMARY OF THE INVENTION

Therefore, in an embodiment of the present invention, a nondestructivetesting method is provided which is capable of analyzing a specificshape of a periodic structure and its internal components.

In a particular embodiment of the present invention, a nondestructivetesting method is provided which is capable of precisely measuring areal shape of a periodic structure, considering a surface layer, such asan oxide layer or a coating layer formed on the periodic structure.

In accordance with one or more embodiments of the present invention,light from a light source is allowed to be incident on a real periodicstructure, so that reflectivity or transmittance of the incident lightor the physical properties related to the reflectivity and transmittanceare measured.

Further, in such embodiment, a virtual periodic structure can be set,and when light is incident on the virtual periodic structure,reflectivity or transmittance of the incident light or the physicalproperties related to the reflectivity and transmittance can becalculated and can be compared with the measured physical properties.

The virtual periodic structure may have periodicity in one dimension,two dimensions or three dimensions. For calculation, the virtualperiodic structure may be horizontally divided into a plurality oflayers. In middle layers of the divided virtual periodic structure,three substances may appear repeatedly in a given period.

The reflectivity, transmittance and the relevant physical properties ofthe virtual periodic structure may be calculated by using refractiveindices of the substances forming the virtual periodic structure. Theactually-measured reflectivity, transmittance and the relevant physicalproperties can be compared with the calculated corresponding physicalproperties and it can be determined whether the formation of the realperiodic structure is identical with that of the virtual periodicstructure.

The reflectivity or transmittance of the virtual periodic structure orboth the reflectivity and the transmittance can be obtained, forexample, from the Fourier series solutions for Maxwell's equationssubject to the so-called Floquet conditions

In accordance with one embodiment of the present invention, the virtualperiodic structure can be divided into a zeroth order structure and aperturbed periodic structure which varies geometrically or physicallyrelative to the zeroth order structure in a perturbation region. Thephysical properties related to the reflectivity and transmittance of thelight incident on the zeroth order structure can be calculated such asby matching Fourier series solutions of Maxwell's equations in alldivided layers, for example. Then, physical properties related to thereflectivity or transmittance of the perturbed virtual periodicstructure or both can be obtained by solving a relevantLippmann-Schwinger equation (See, e.g., Eq. 9 below), with Green'sfunction (See, e.g., Eq. 8_below) of the zeroth order structure.

Accordingly, in accordance with one or more embodiments of the presentinvention, the formation and internal components of the periodicstructure can be non-destructively tested.

In accordance with a particular embodiment of the present invention,microscopic change relative to an original periodic structure having anative oxide layer thereon or intentionally formed surface coating layerthereon can be precisely tested. The development of the semiconductorindustry or other nano technology can benefit from such embodiment ofthe invention.

BRIEF DESCRIPTION OF THE DRAWINGS

Features and advantages of embodiments of the present invention will beapparent in accordance with the detailed description of preferredembodiments thereof with reference to the attached drawings in which:

FIG. 1 is a flow chart schematically illustrating a testing methodaccording to an embodiment of the present invention;

FIGS. 2 and 3 are block diagrams schematically illustrating testingdevices of a periodic structure;

FIG. 4 is a perspective view illustrating an example of a virtualperiodic structure;

FIG. 5 is a sectional view of the virtual periodic structure of FIG. 4being divided into a plurality of layers;

FIG. 6 is a perspective view of the geometrical formation of a virtualperiodic structure according to an embodiment of the present invention;

FIG. 7 is a sectional view of the virtual periodic structure of FIG. 6being divided into a plurality of layers;

FIG. 8 is a sectional view of the geometrical formation of a virtualperiodic structure according to an embodiment of the present invention;

FIGS. 9 and 10 are graphs respectively illustrating results ofcalculating the reflectivity of a virtual periodic structure consideringan oxide layer and calculating the reflectivity of a virtual periodicstructure not considering the oxide layer in a TE mode and a TM mode;and

FIGS. 11 and 12 are graphs illustrating results of extracting thephysical properties measured by ellipsometry from the calculationresults of FIGS. 9 and 10.

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

Particular embodiments of the present invention will now be describedmore fully hereinafter with reference to the accompanying drawings.

FIG. 1 is a flow chart schematically illustrating a testing methodaccording to an embodiment of the present invention. By the testingmethod illustrated in FIG. 1, the formation and internal component of aperiodic structure can be non-destructively tested. Specifically, shape,dimensions or both shape and dimensions of the periodic structure can beprecisely determined, even when the periodic structure includes a nativeoxide layer formed on the periodic structure or an intentionally formedsurface coating layer.

Measurement of Periodic Structure

FIGS. 2 and 3 are block diagrams schematically illustrating testingdevices for a periodic structure 200. Referring to FIG. 2, the testingdevices comprise a light source 100, a detector 110 and a processor 120.When a periodic structure as a test object is positioned on a substrate130, the light source 100 irradiates light having a specific wavelengthor various wavelengths to the periodic structure 200. The light incidenton the periodic structure 200 is partially transmitted and partiallyreflected. The reflected light is detected in the detector 110, and thereflectivity of a reflected wave measured in the detector 110 iscalculated in the processor 120. The transmitted light is also detectedin the detector 110, and the transmittance of a transmitted wavemeasured in the detector 110 is calculated in the processor 120.

The testing device may further comprise a polarizer 140 as illustratedin FIG. 3. In this case, the light generated from the light source 100is polarized as the light of a TE mode or TM mode through the polarizer140, to be incident on the periodic structure 200.

When the light is incident on the periodic structure 200, the incidentlight is divided into the reflected light and the transmitted light. Inan embodiment of the present invention, reflectivity and transmittancein the most basic two polarization states in the reflection andtransmission of the light, that is, the TE mode and the TM mode, arecalculated to perform the nondestructive test of the periodic structure.

For example, the physical properties related to the reflectivity andtransmittance which are measured by allowing the light to be incident onthe periodic structure may be understood as the combination of thephysical properties, which are related to an amplitude or phase of areflected wave and a transmitted wave to an incident wave of a TE modeelectric field, and the physical properties, which are related to anamplitude or phase of a reflected wave and a transmitted wave to anincident wave of a TM mode magnetic field.

Virtual Periodic Structure

FIG. 4 is a perspective view illustrating the geometrical formation of avirtual periodic structure 200 a, and FIG. 5 illustrates that a sectionof the virtual periodic structure 200 a of FIG. 4 is divided into aplurality of layers.

In a case where the virtual periodic structure 200 a is a semiconductordevice having a shape of, for example, one-dimensional, two-dimensionalor three-dimensional periodic formation, the virtual periodic structure200 a has the formation in which two substances of substance parts (fromn₁ to n_(L)) formed of a substance, such as silicon and the like, ineach layer and substance parts (from ñ₁ to ñ₁) of an incident part, suchas an air layer and the like, are horizontally periodic.

However, in the real environments of the semiconductor process and thelike, since it is impossible to fabricate the periodic structure in aperfect vacuum state, an oxide layer is formed on the surface of thereal periodic structure by contact with air or water. Further, in theprocessing steps, since an intentional coating layer is formed on thesurface of the periodic structure or a roughness layer is present on thesurface of the periodic structure, the periodic structure of FIG. 4 hasa limit in presuming the real geometrical shape thereof.

FIG. 6 illustrates the geometrical aspects, e.g., shape of a virtualperiodic structure 200 b according to an embodiment of the presentinvention. On the surface of the virtual periodic structure 200 b, asurface layer 210, such as an oxide layer and the like, is formed. Whena section of the virtual periodic structure 200 b is divided into aplurality of layers, the virtual periodic structure 200 b has theformation in which at least three substances are periodically repeated,as illustrated in FIG. 7. The virtual periodic structure 200 b includesa ridge region formed at both sides of a third substance correspondingto a groove region. This ridge region is formed of a first substanceforming a center part, and a second substance. The second substanceincludes the surface layer formed on the outer surface of the centerpart. In FIG. 7, n_(l) (l=2, . . . , L), ñ_(l) (l=1, . . . , L−1), and n_(l) (l=1, . . . , L) are respectively representative of the ridgeregion (the first substance), the groove region (the third substance)and the surface layer region (the second substance). The ridge regionmay be formed of a single substance or two or more different substancesbeing vertically or horizontally positioned. Accordingly, the firstsubstance of the ridge region may be formed of a plurality ofsubstances.

The surface layer region (the second substance) may be the oxide layeror the coating layer, or it may be the roughness layer of the periodicstructure surface as the case may be. A third substance occupying thegroove region may be in gaseous, liquid or solid state or a combinationthereof.

For example, when the virtual periodic structure 200 b is asemiconductor device, the plurality of layers (1 to L) except for thehighest layer (layer 1) can be formed of the first substance being asemiconductor such as silicon and the like and the highest layer (layer1) being a the second substance such as an oxide layer or a coatinglayer. A third substance being an air layer or other gas, liquid orsolid state substance can be disposed in grooves between the layers 1 toL. The layers 1 to L and the third substance can be horizontally andperiodically repeated within the periodic structure.

When the virtual periodic structure 200 b is set by considering thesurface layer 210, such as the oxide layer, coating layer or roughnesslayer of the surface, the reflectivity and transmittance of the virtualperiodic structure is calculated to more closely approximate those ofthe real periodic structure. Consequently, the shape and components ofthe real periodic structure can be accurately measured. Specifically, itis possible to compare and analyze the geometrical shape and theinternal component of the periodic structure, including the thicknessesof thin film structures present within the periodic structure.

Calculation of Physical Properties of Virtual Periodic Structure

The reflectivity and transmittance of light by a virtual periodicstructure can be calculated by expanding refractive indices in Fourierseries, and electric and magnetic fields in Floquet series. To accountfor vertical variation of the refractive index in the verticaldirection, the direction aligned with the z axis (FIG. 4), the virtualperiodic structure 200 b is divided into the plurality of layers (1 toL) illustrated in FIG. 7, the refractive of index of each layer isassumed to be constant over the vertical height of the layer. Moreprecise solutions can be obtained by increasing the number of the layersdivided from the virtual periodic structure 200 b and the number ofterms of the Fourier series expansion.

The complex expansion coefficients of electric and magnetic waves inregion I and region II gives the amplitude and phase of the reflectedand transmitted waves, respectively. The expansion coefficients ofdielectric function in each layer can be expressed in terms of complexrefractive index of each of the first and second substances forming theridge region and the third substance forming the groove region, a ratiof, of the region occupied by the first and second substances to a periodΛ, a ratio 2δ₁/Λ of the region occupied only by the second substance,and a parameter t_(l) indicating how far the center of the layer l isfrom the center of layer 1 in the x direction, which will be describedlayer.

Further, in an embodiment of the present invention, to calculate thephysical properties of the virtual periodic structure, a shape of thevirtual periodic structure is assumed as a zeroth order periodicstructure. Then, a perturbed periodic structure is obtained by applyinga geometrical or physical change to the zeroth order periodic structurein a perturbation region.

To obtain the reflectivity or transmittance of the perturbed periodicstructure, the Lippmann-Schwinger equation (Eq. 9 below) is solvednumerically using the numerically calculated Green's function,reflectivity, and transmittance of the zeroth order structure as inputvalues.

As described above, the reflectivity and transmittance of the zerothorder periodic structure can be calculated by solving Maxwell's Eqs.with Floquet's condition. For this purpose, the virtual periodicstructure 200 b is divided into the plurality of layers (1 to L) in therectangular sectional shape, and subsequently, the dielectric functionin each layer is expanded in Fourier series. When light is incident onthe virtual periodic structure 200 b, the reflected and transmittedwaves in each layer are also expanded in Floquet series withcoefficients to be determined by matching conditions of the electric andmagnetic fields.

From the expansion coefficients thus determined, the reflectivity(ρ_(TE) ^([0]) and ρ_(TM) ^([0])) and the transmittance (τ_(TE) ^([0])and τ_(TM) ^([0])) of the zeroth order periodic structure arecalculated.

The Green's function of the zeroth order structure in the perturbationregion can be obtained by applying the method of calculating thetransmittance or reflectivity of the zeroth order periodic structure.

The reflectivity (ρ_(TE) and ρ_(TM)) and the transmittance (τ_(TE) andτ_(TM)) of the perturbed periodic structure can be calculated from theexpansion coefficients of a coupled wave of the TE mode electric fieldand the TM mode magnetic field in each layer of perturbed structureobtained by solving the Lippmann-Schwinger equation.

In accordance with an embodiment of the present invention, the periodicstructure between an incident medium (region I) and a substrate (regionII) is divided into L layers, including some of or the whole of layersbeing formed of a uniform substance. The divided grating structures areclassified into the following cases: the case where each layer in alllayers of the whole perturbation region is formed of a uniform substance(case 1); (ii) the case where all layers of the perturbation region areformed of a single uniform substance (case 2); (iii) the case where theperturbation region is placed as the air layer on the substrate (case3); and (iv) the case where each layer is not formed of a uniformsubstance but all layers are formed of one given layer being repeated(case 4).

In accordance with an embodiment of the present invention, there areprovided the results of specific calculation of each case 1, case 2,case 3 and case 4, specifically, considering the periodic structure withthe surface layer, which will be described later.

Comparison of the Measured Value with the Calculated Value

The actually-measured reflectivity or transmittance of the periodicstructure are compared with the calculated reflectivity or transmittanceof the virtual periodic structure. Alternatively, the actually-measuredreflectivity or transmittance of the periodic structure are comparedwith the calculated reflectivity or transmittance of the virtualperiodic structure. In another alternative, properties related to theactual reflectivity or transmittance of the periodic structure arecompared with properties related to the calculated reflectivity ortransmittance of the virtual periodic structure. When these values arethe same within a given error range, the real periodic structure, asdetermined from measurements of reflectivity or transmittance, can bedetermined as being identical with that of the virtual periodicstructure 200 b determined by calculating values of reflectivity ortransmittance.

When comparing the measured reflectivity and transmittance with thecalculated reflectivity and transmittance, an additional device, suchas, for example, a computer may be used for comparing the measured valuewith the calculated value. By this method, geometry, e.g., shape anddimensional aspects of the real periodic structure 200 can be preciselyand efficiently determined.

In the measurement experiment of the periodic structure, diversephysical properties related to reflectivity ρ and transmittance τ aremeasured. The relevant presumed value related to the measured physicalproperties can be obtained by the reflectivity ρ and transmittance τbeing calculated based on the mathematical methods.

Embodiment

As an embodiment of the present invention, a virtual periodic structurewhere a surface layer is formed is set as illustrated in FIG. 8. Thisvirtual periodic structure has a repeated period of 60 nm. The thicknessof a ridge region formed of silicon crystals is 45 nm, an oxide isformed on the surface of the ridge region to the thickness of 5 nm, andthe slope of a ridge is set as 81°.

The slope and period of the virtual periodic structure and the width andheight of the ridge may be varied as random values. Further, thegeometrical shape thereof may have a random shape.

The physical properties of the virtual periodic structure are calculatedby using the method according to an embodiment of the present invention,to be compared with the measured values of a real periodic structure.

FIGS. 9 and 10 are graphs respectively illustrating results of thecalculated reflectivity of the virtual periodic structure consideringthe oxide layer as illustrated in FIG. 8 and the calculated reflectivityof the virtual periodic structure not considering the oxide layer in aTE mode and a TM mode. That is, the reflectivity of an incidentwavelength in each of the TE mode and the TM mode is calculated,assuming a case where the oxide layer is considered, i.e., the virtualperiodic structure is formed of silicon coated with the oxide layer, andanother case where the oxide layer is not considered, i.e., the virtualperiodic structure is formed of pure silicon.

Based on the results as illustrated, the calculated value of thereflectivity of the pure silicon is considerably different from that ofthe reflectivity of the silicon coated with the oxide layer. That is,the extent of approximation of the virtual periodic structure to thereal periodic structure significantly varies, depending on whether ornot we assume the presence of the surface layer, such as the oxide layerand the like. In the case of measuring the formation of the realperiodic structure in which the surface layer, such as the oxide layerand the like, is formed, if the virtual periodic structure is assumed tobe formed of the pure silicon only, the calculated values of the virtualperiodic structure are never identical with the measured values of thereal periodic structure. When the real periodic structure to be measuredincludes the surface layer, such as the oxide layer and the like, theextent of identity with the real periodic structure obviously increasesonly when the virtual periodic structure is assumed to be formed ofsilicon coated with the oxide layer.

FIGS. 11 and 12 are graphs illustrating results of the extractedphysical properties measured by ellipsometry from the calculationresults of FIGS. 9 and 10. The graphs show that diverse physicalproperties can be obtained from the results of calculation using themethod in accordance with an embodiment of the present invention.Specifically, the graphs illustrate the results of the calculatedamplitude and phase of the reflectivity ratio between the TM mode ρ_(M)and the TE mode ρ_(E) as the physical properties Ψ and Δ measured bygeneral ellipsometry. The analysis results of the graphs show that thereis a big difference between the virtual periodic structure formed of thepure silicon and the virtual periodic structure formed of the siliconcoated with the oxide layer. As described above, upon comparing with theformation of the real periodic structure by extracting the diversephysical properties from the reflectivity or transmittance, if theformation of the virtual periodic structure is assumed considering thepresence of the oxide layer, the more accurate formation of the realperiodic structure is presumed.

The method of calculating the reflectivity and transmittance of thevirtual periodic structure 200 b illustrated in FIGS. 6 and 7 accordingto other embodiments of the present invention will be described below:

I. Lippmann-Schwinger Equation: TE Mode

Due to the periodicity of the periodic structure in x direction, the TEmode solutions to Maxwell's equation can be written as

$\begin{matrix}{{E_{y}\left( {x,z} \right)} = {\sum\limits_{n = {- \infty}}^{\infty}\; {{S_{yn}(z)}^{{k}_{xn}x}}}} & {{Eq}.\mspace{14mu} (1)} \\{{H_{x}\left( {x,z} \right)} = {{i\left( \frac{\in_{0}}{\mu_{0}} \right)}^{1/2}{\sum\limits_{n = {- \infty}}^{\infty}\; {{U_{xn}(z)}^{{k}_{xn}x}}}}} & {{Eq}.\mspace{14mu} (2)}\end{matrix}$

where i=√{right arrow over (−1)}; k_(xn)=k₀[n₁ sin θ−n(λ₀/Λ)]; k₀=2π/λ₀(with λ₀ being the wavelength of an incident wave in vacuum, n₁ beingthe refractive index of an incident region, and θ being the incidentangle); ε₀ is the permittivity of the vacuum, and μ₀ is the permeabilityof the vacuum.

The dielectric function of the periodic structure with a one-dimensionalperiod being equal to Λ in x direction can be expanded in Fourierseries:

$\begin{matrix}{{ɛ\left( {x,z} \right)} = {\sum\limits_{h = {- \infty}}^{\infty}\; {{ɛ_{h}(z)}^{{2}\; \pi \; {{hx}/\Lambda}}}}} & {{Eq}.\mspace{14mu} (3)}\end{matrix}$

From Maxwell's equations, S_(yn)(z) in Eq. (1) and U_(xn)(z) in Eq. (2)are combined as follows:

$\begin{matrix}{\begin{bmatrix}{\left( {1/k_{0}} \right){{\partial{S_{y}(z)}}/{\partial z}}} \\{\left( {1/k_{0}} \right){{\partial{U_{x}(z)}}/{\partial z}}}\end{bmatrix} = {\begin{bmatrix}0 & I \\{A(z)} & 0\end{bmatrix}\begin{bmatrix}{S_{y}(z)} \\{U_{x}(z)}\end{bmatrix}}} & {{Eq}.\mspace{14mu} (4)}\end{matrix}$

where I is a unit matrix of dimension 2N+1 and

A(z)≡K _(x) ² −E(z)  Eq. (5)

where K_(x) is a diagonal matrix with the (n,n) element being equal tok_(xn)/k₀ and E(z) is a matrix formed by the permittivity harmoniccomponents, with the (n,p) element being equal toE_(np)(z)=ε_((n−p))(z).

Below, coupled-wave-basis indices are truncated so that they run from −Nto N.

Combining two first order differential equations in Eq. (4), we obtain

$\begin{matrix}{{\left\lbrack {{\frac{1}{k_{0}^{2}}\frac{\partial^{2}}{\partial z^{2}}} - {A(z)}} \right\rbrack {S_{y}(z)}} = 0} & {{Eq}.\mspace{14mu} (6)}\end{matrix}$

With the superscript [0], indicating the zeroth order structure, weobtain a similar Eq.:

$\begin{matrix}{{\left\lbrack {{\frac{1}{k_{0}^{2}}\frac{\partial^{2}}{\partial z^{2}}} - {A^{\lbrack 0\rbrack}(z)}} \right\rbrack {S_{y}^{\lbrack 0\rbrack}(z)}} = 0} & {{Eq}.\mspace{14mu} (7)}\end{matrix}$

A corresponding Green's function can be formulated as follows:

$\begin{matrix}{{\left\lbrack {{\frac{1}{k_{0}^{2}}\frac{\partial^{2}}{\partial z^{2}}} - {A^{\lbrack 0\rbrack}(z)}} \right\rbrack {G\left( {z,z^{\prime}} \right)}} = {{- {\delta \left( {z - z^{\prime}} \right)}}I}} & {{Eq}.\mspace{14mu} (8)}\end{matrix}$

From Eqs. (6), (7) and (8) are, we obtain the Lippman-Schwingerequation:

S _(y)(z)−S _(y) ^([0])(z)=∫dz′G(z,z′)  Eq. (9)

where Y(z)≡E(z)−E(z)−E^([0])(z). In Eq. (9), since Y(z)=O outside theperturbation zone, the integration region is restricted within theperturbation region.

II. Calculation of S_(y) ^([0])(z)

In this section, we omit the superscript [0] used in the previoussection to indicate the relevant quantities in the zeroth order periodicstructure.

A. Division of Periodic Structure Region into Layers

We approximate the periodic structure, which includes an oxide layer, acoating layer or a surface layer, to a stack of parallel rectangularlayers with a common period Λ, whereby the z dependency of thedielectric function ε(x,z) is shifted to the index l indicating thelayer. Then, the dielectric function of the given layer l is indicatedas ε^((l))(x) and can be expressed as follows:

$\begin{matrix}{{ɛ^{(l)}(x)} = {\sum\limits_{h}{ɛ_{h}^{(l)}^{\; 2\; \pi \; {{hx}/\Lambda}}}}} & {{Eq}.\mspace{14mu} (10)}\end{matrix}$

The expansion coefficient ε_(h) ^((l)) in Eq. (10) is given in eachlayer as follows:

$\begin{matrix}{l = 1} & \; \\{ɛ_{0}^{(1)} = {{{\overset{\_}{n}}_{1}^{2}f_{1}} + {{\overset{\sim}{n}}_{1}^{2}\left( {1 - f_{1}} \right)}}} & {{Eq}.\mspace{14mu} (11)} \\{ɛ_{h}^{(1)} = {{\left( {{\overset{\_}{n}}_{1}^{2} - {\overset{\sim}{n}}_{1}^{2}} \right)\frac{\sin \left( {\pi \; {hf}_{1}} \right)}{\pi \; h}\mspace{14mu} {from}\mspace{14mu} l} = {{2\mspace{14mu} {to}\mspace{14mu} L} - 1}}} & {{Eq}.\mspace{14mu} (12)} \\{ɛ_{0}^{(l)} = {\left\lbrack {{\left( {n_{l}^{2} - {\overset{\sim}{n}}_{l}^{2}} \right)f_{l}} + {\left( {{\overset{\_}{n}}_{l}^{2} - n_{l}^{2}} \right)\frac{2\delta_{l}}{\Lambda}} + {\overset{\sim}{n}}_{l}^{2}} \right\rbrack ^{{- }\; 2\; \pi \; {L_{i}/\Lambda}}}} & {{Eq}.\mspace{14mu} (13)} \\{\begin{matrix}{ɛ_{h}^{(l)} = {\frac{1}{\pi \; h}\left\lbrack {{\left( {{\overset{}{n}}_{l}^{2} - {\overset{\sim}{n}}_{l}^{2}} \right){\sin \left( {\pi \; {hf}_{l}} \right)}} +} \right.}} \\\left. {\left( {n_{l}^{2} - {\overset{\_}{n}}_{l}^{2}} \right){\sin \left( {{\pi \; {hf}_{i}} - \frac{2\; \pi \; h\; \delta_{l}}{\Lambda}} \right)}} \right\rbrack\end{matrix}^{{- }\; 2\; \pi \; {t_{l}/\Lambda}}} & {{Eq}.\mspace{14mu} (14)} \\\begin{matrix}{ɛ_{0}^{(l)} = \left\lbrack {{n_{L}^{2}\left( {f_{L - 1} - \frac{2\delta_{L - 1}}{\Lambda}} \right)} +} \right.} \\{\left. {{\overset{\_}{n}}_{L}^{2}\left( {1 - f_{L - 1} + \frac{2\delta_{L - 1}}{\Lambda}} \right)} \right\rbrack ^{{- }\; 2\; \pi \; {t_{L - 1}/\Lambda}}}\end{matrix} & {{Eq}.\mspace{14mu} (15)} \\{ɛ_{h}^{(l)} = {\frac{1}{\pi \; h}\left( {n_{L}^{2} - {\overset{\_}{n}}_{L}^{2}} \right){\sin \left\lbrack {\pi \; {h\left( {f_{L - 1} - \frac{2\delta_{L - 1}}{\Lambda}} \right)}} \right\rbrack}^{{- {2}}\; \pi \; {t_{L - 1}/\Lambda}}}} & {{Eq}.\mspace{14mu} (16)}\end{matrix}$

where n_(l) (l=2, . . . , L), ñ_(l) (l=1, . . . , L−1), and n _(l) (l=1,. . . , L) are respectively complex refractive index of the ridgeregion, groove region and surface layer region in the layer l, f^(l) isthe fraction of the period Λ occupied by the ridge region (the first andsecond substances), 2δ_(l)/Λ is the fraction of the period occupied bythe oxide layer region (the second substance) only, and t_(l) (t₁=t₂=O)is the center-shifting parameter of the layer l relative to the centerof layer 1 in x direction.

In the division of the periodic structure into layers, some of or wholeof the layers may be formed of a uniform substance. As described above,the multilayer Ga structure indicates the case where each layer in alllayers is formed of a uniform substance, and the single layer Gastructure indicates the case where all layers are formed of a singleuniform substance. Further, when each layer is not formed of a uniformsubstance but all layers are formed of one given layer which isrepeated, it is defined as a single layer G1 structure. Specifically,region II of FIG. 6 may be an air layer.

The solutions of the electric field and magnetic field in regions I andII can be represented by using Floquet's conditions as follows:

$\begin{matrix}{{E_{y}^{(0)}\left( {x,z} \right)} = {\sum\limits_{n = {- N}}^{N}{{S_{yn}^{(0)}(z)}^{{k}_{xn}x}}}} & {{Eq}.\mspace{14mu} (17)} \\{{H_{x}^{(0)}\left( {x,z} \right)} = {{\left( \frac{\varepsilon_{0}}{\mu_{0}} \right)}^{1/2}{\sum\limits_{n = {- N}}^{N}{{U_{xn}^{(0)}(z)}^{{k}_{xn}x}}}}} & {{Eq}.\mspace{14mu} (18)} \\{{E_{y}^{({L + 1})}\left( {x,z} \right)} = {\sum\limits_{n = {- N}}^{N}{{S_{yn}^{({L + 1})}(z)}^{{k}_{xn}x}}}} & {{Eq}.\mspace{14mu} (19)} \\{{H_{x}^{({L + 1})}\left( {x,z} \right)} = {\; \left( \frac{\varepsilon_{0}}{\mu_{0}} \right)^{1/2}{\sum\limits_{n = {- N}}^{N}{{U_{xn}^{({L + 1})}(z)}^{{k}_{xn}x}}}}} & {{Eq}.\mspace{14mu} (20)} \\{{S_{yn}^{(0)}(z)} = {{\delta_{0n}^{{k}_{L,{zn}}z}} + {\rho_{n}^{{- {k}_{L,{zn}}}z}}}} & {{Eq}.\mspace{14mu} (21)} \\{{U_{xn}^{(0)}(z)} = {{{\delta_{0n}\left( {{k}_{L,{zn}}/k_{0}} \right)}^{{k}_{L,{zn}}z}} - {{\rho_{n}\left( {{k}_{L,{zn}}/k_{0}} \right)}^{{- {k}_{L,{zn}}}z}}}} & {{Eq}.\mspace{14mu} (22)} \\{{S_{yn}^{({L + 1})}(z)} = {\tau_{n}^{{k}_{11,{zn}}{({z - Z_{L}})}}}} & {{Eq}.\mspace{14mu} (23)} \\{{U_{xn}^{({L + 1})}(z)} = {{\tau_{n}\left( {{k}_{11,{zn}}/k_{0}} \right)}^{{k}_{11,{zn}}z}}} & {{Eq}.\mspace{14mu} (24)} \\{k_{,{zn}} = \left\{ \begin{matrix}{k_{0}\left\lbrack {n_{}^{2} - \left( {k_{xn}/k_{0}} \right)^{2}} \right\rbrack}^{1/2} & {{{for}\mspace{14mu} k_{0}n_{}} > k_{xn}} \\{{k}_{0}\left\lbrack {\left( {k_{xm}/k_{0}} \right)^{2} - n_{}^{2}} \right\rbrack}^{1/2} & {{{{{for}\mspace{14mu} k_{xn}} > {k_{0}n_{}}};{ = I}},{II}}\end{matrix} \right.} & {{Eq}.\mspace{14mu} (25)}\end{matrix}$

where n_(I) and n_(II) are the complex refractive index of the incidentmedium and the substrate. The first term on the right-hand side of Eq.(21) represents the incident wave.

The electric field and the magnetic field in the periodic structureregions (l=1, . . . , L) can be expanded in terms of coupled-wave bases:

$\begin{matrix}{{E_{y}^{(l)}\left( {x,z} \right)} = {\sum\limits_{n = {- N}}^{N}{{S_{yn}^{(l)}(z)}^{{k}_{xn}x}}}} & {{Eq}.\mspace{14mu} (26)} \\{{H_{x}^{(l)}\left( {x,z} \right)} = {{\left( \frac{\varepsilon_{0}}{\mu_{0}} \right)}^{1/2}{\sum\limits_{n = {- N}}^{N}{{U_{xn}^{(l)}(z)}^{{k}_{xn}x}}}}} & {{Eq}.\mspace{14mu} (27)}\end{matrix}$

When Eq. (7) is applied to the layer l, we obtain an Eq. for S_(yn)^((l)) (n=−N, . . . , N) in Eq. (26) as follows:

(1/k ₀ ²)∂² S _(y) ^((l))(z)/∂z ² =A _(l) S _(y) ^((l))(z)  Eq. (28)

where

A _(l) ≡K _(x) ² −E _(l)  Eq. (29)

where the (n,p) element of E_(l) is equal to

Solutions S_(y) ^((l))(z) and U_(x) ^((l))(z) to Eq. (28) are given asfollows:

S _(y) ^((l))(z)=W _(l) [e ^(−k) ⁰ ^(Q) ^(l) ^((z−z) ^(l−1) ⁾ f _(l) +e^(k) ⁰ ^(Q) ^(l) ^((z−z) ^(l−1) ⁾ g _(l)]  Eq. (30)

U _(x) ^((l))(z)=V _(l) [−e ^(−k) ⁰ ^(Q) ^(l) ^((z−z) ^(l−1) ⁾ f _(l) +e^(k) ⁰ ^(Q) ^(l) ^((z−z) ^(l−1) ⁾ g _(l)]  Eq. (31)

where W_(l) is a square matrix constructed from (2N+1) eigenvectors ofA_(l), Q_(l) is a diagonal matrix with the element being the positivesquare root of (2N+1) eigenvalue of A_(l), and V_(l)=W_(l)Q_(l); f_(l)and g_(l) are column vectors to be determined from boundary conditions;z_(l−1) is the z coordinate of the interface between the layer l−1 andthe layer l.

Eqs. (21), (22), (23) and (24), respectively, can be written in thefollowing matrix form:

S _(y) ⁽⁰⁾(z)=W ₀ [e ^(−k) ⁰ ^(Q) ⁰ ^(z) f ₀ +e ^(k) ⁰ ^(Q) ⁰ ^(z) g₀],  Eq. (32)

U _(x) ⁽⁰⁾(z)=V ₀ [−e ^(−k) ⁰ ^(Q) ⁰ ^(z) f ₀ +e ^(k) ⁰ ^(Q) ⁰ ^(z) g₀],  Eq. (33)

S _(y) ^((L+1))(z)=W _(L+1) [e ^(−k) ⁰ ^(Q) ^(L+1) ^((z−z) ^(L) ⁾ f_(L+1) +e ^(k) ⁰ ^(Q) ^(L+1) ^((z−z) ^(L) ⁾ g _(L+1)],  Eq. (34)

U _(x) ^((L+1))(z)=V _(L+1) [−e ^(−k) ⁰ ^(Q) ^(L+1) ^((z−z) ^(L) ⁾ f_(L+1) +e ^(k) ⁰ ^(Q) ^(L+1) ^((z−z) ^(L) ⁾ g _(L+1)]  Eq. (35)

where Q₀ and Q_(L+1) are diagonal matrices with the n-th diagonalelements being −ik_(I,zn)/k₀ and −ik_(II,zn)/k₀, respectively; all off₀=(0, . . . , 0, 1, 0, . . . , 0)^(T), g₀=(ρ_(−N), . . . ρ⁻¹, ρ₀, ρ₁, .. . ρ_(N))^(T), f_(L+1)=(τ_(−N), . . . , τ⁻¹, τ₀, τ₁, . . . ,τ_(N))^(T), g_(L+1)=(0, . . . , 0)^(T) are (2−N+1)-component columnvectors; W₀=I, V₀=W₀Q₀, W_(L+1)=I and V_(l+1)=W_(L+1)Q_(L+1).are

B. Calculation of f_(l) and g_(l)

By applying the boundary conditions which the electric field and themagnetic field satisfy at z=z_(l) (l=0, . . . , L), we have

$\begin{matrix}{{\begin{bmatrix}f_{l} \\g_{l}\end{bmatrix} = {{\begin{bmatrix}X_{l}^{- 1} & 0 \\0 & X_{l}\end{bmatrix}\begin{bmatrix}M_{l}^{+} & M_{l}^{-} \\M_{l}^{-} & M_{l}^{+}\end{bmatrix}}\begin{bmatrix}f_{l + 1} \\g_{l + 1}\end{bmatrix}}},{where}} & {{Eq}.\mspace{14mu} (36)} \\{X_{l} \equiv \left\{ {\begin{matrix}^{{- k_{0}}Q_{i}d_{l}} & {{{{for}\mspace{14mu} l} = 1},\ldots \mspace{14mu},L} \\I & {{{for}\mspace{14mu} l} = 0}\end{matrix},} \right.} & {{Eq}.\mspace{14mu} (37)} \\{M_{l}^{\pm} \equiv {\frac{1}{2}\left\lbrack {{W_{l}^{- 1}W_{l + 1}} \pm {V_{l}^{- 1}V_{1 + 1}}} \right\rbrack}} & {{Eq}.\mspace{14mu} (38)}\end{matrix}$

Because A_(l) is symmetric matrix, eigenvectors of A_(l) are orthogonal.Thus, when the normalized-eigenvector matrix W_(l) is used, W_(l) ^(T)for W_(l) ⁻¹ may be used to enhance the efficiency of numericalcalculation.

The recursion relation of Eq. (36) will eventually yield the followingform:

$\begin{matrix}{\begin{bmatrix}f_{l} \\g_{l}\end{bmatrix} = {\begin{bmatrix}a_{l} \\b_{l}\end{bmatrix}f_{L + 1}}} & {{Eq}.\mspace{14mu} (39)}\end{matrix}$

Eliminating f_(L+1) from Eq. (39), we obtain:

g_(l)=b_(l)a_(l) ⁻¹f_(l)≡R_(l)f_(l)  Eq. (40)

Substituting these g_(l) and g_(l+1)=R_(l) f _(l+1) into Eq. (36), weobtain the recursion relation for R_(l) (l=0, 1, . . . , L):

R _(l) =X _(l)(M_(l) ⁻ +M _(l) ⁺ R _(l+1))(M _(l) ⁺ +M _(l) ⁻ R_(l+1))⁻¹ X _(l).  Eq. (41)

By applying this recursion relation repeatedly starting from l=L withthe initial setting value R_(L+1)=O obtain all values of R_(l).

Now that we have an algorithm to determine R_(l) (l=0, . . . , L) in Eq.(40), the remaining problem in determining S^((l)) s completely is tocalculate f_(l)s. Since f₀ is known, f_(l) can be calculated by applyingthe transfer matrix method from the incident medium to layer l.

To this end, the recursion relation of Eq. (36) is inverted to obtain

$\begin{matrix}{f_{l + 1} = {\left( {{N_{l}^{+}X_{l}} + {N_{l}^{-}X_{l}R_{l}}} \right)f_{l}\mspace{14mu} {where}}} & {{Eq}.\mspace{14mu} (42)} \\{N_{l}^{\pm} = {\frac{1}{2}\left\lbrack {{W_{l + 1}^{- 1}W_{l}} \pm {V_{l + 1}^{- 1}V_{l}}} \right\rbrack}} & {{Eq}.\mspace{14mu} (43)}\end{matrix}$

By repeated applications of Eq. (42) with the initial setting f₀=(0, . .. , 0, 1, 0, . . . , 0)^(T), all f_(l)s are determined.

III. Calculation of G(z,z′): TE Mode

If z locates inside the layer l, Eq. (8) becomes:

$\begin{matrix}{{\left\lbrack {{\frac{1}{k_{0}^{2}}\frac{\partial^{2}}{\partial z^{2}}} - A_{l}} \right\rbrack {G\left( {z,z^{\prime}} \right)}} = {{- {\delta \left( {z - z^{\prime}} \right)}}I}} & {{Eq}.\mspace{14mu} (44)}\end{matrix}$

In this section, we omit, as in section II, the superscript [0] used toindicate the relevant quantities in the zeroth order periodic structure.For the given layer l, A^(l) is not a function of z any more. Thesolutions to Eq. (44) are of the following forms:

$\begin{matrix}{G_{{ll}^{\prime}{({z,z^{\prime}})}} = \left\{ \begin{matrix}\begin{matrix}\begin{matrix}{W_{l}\left\lbrack {{^{{- k_{0}}{Q_{l}{({z - z_{l - 1}})}}}f_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)} +} \right.} \\\left. {^{k_{0}{Q_{l}{({z - z_{l - 1}})}}}{g_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)}} \right\rbrack\end{matrix} \\\left( {{{for}\mspace{14mu} l} < l^{\prime}} \right)\end{matrix} & \; \\\begin{matrix}\begin{matrix}{W_{l}\left\lbrack {{^{{- k_{0}}{Q_{l}{({z - z_{l - 1}})}}}f_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)} +} \right.} \\\left. {^{k_{0}{Q_{l}{({z - z_{l - 1}})}}}{g_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)}} \right\rbrack\end{matrix} \\\left( {{{{for}\mspace{14mu} l} = l^{\prime}};{z < z^{\prime}}} \right)\end{matrix} & \; \\\begin{matrix}\begin{matrix}{W_{l}\left\lbrack {{^{{- k_{0}}{Q_{l}{({z - z_{l - 1}})}}}f_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)} +} \right.} \\\left. {^{k_{0}{Q_{l}{({z - z_{l - 1}})}}}{g_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)}} \right\rbrack\end{matrix} \\\left( {{{{for}\mspace{14mu} l} = l^{\prime}};{z > z^{\prime}}} \right)\end{matrix} & \; \\\begin{matrix}\begin{matrix}{W_{l}\left\lbrack {{^{{- k_{0}}{Q_{l}{({z - z_{l - 1}})}}}f_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)} +} \right.} \\\left. {^{k_{0}{Q_{l}{({z - z_{l - 1}})}}}{g_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)}} \right\rbrack\end{matrix} \\\left( {{{for}\mspace{14mu} l} > l^{\prime}} \right)\end{matrix} & \;\end{matrix} \right.} & {{Eq}.\mspace{14mu} (45)}\end{matrix}$

where W_(l) and Q_(L) are square matrices constructed from eigenvectorsand eigenvalues of A_(l) as in section II; and f_(ll′) ^(±)(z′) andg_(ll′) ^(±)(z′), which are constant (in z) square matrices to bedetermined by matching the boundary conditions; and the pairedsubscripts u′ in G_(ll′)(z,z′), f_(ll′) ^(±)(z′), and g_(ll′) ^(±)(z′),are to denote that a field point z locates inside the layer l and asource point z′ locates inside the layer l′. When an index l=O isassigned to region I and an index l=L+1 is assigned to region II, thesolutions in these regions can be written as follows:

$\begin{matrix}{G_{{ll}^{\prime}{({z,z^{\prime}})}}\left\{ \begin{matrix}\begin{matrix}\begin{matrix}{W_{0}\left\lbrack {{^{{- k_{0}}Q_{0z}}f_{0l^{\prime}}^{-}\left( z^{\prime} \right)} +} \right.} \\\left. {^{k_{0}Q_{0z}}{g_{0l^{\prime}}^{-}\left( z^{\prime} \right)}} \right\rbrack\end{matrix} \\\left( {{{for}\mspace{14mu} l} = {0 < l^{\prime}}} \right)\end{matrix} & \; \\\begin{matrix}{W_{L + 1}\left\lbrack {{^{{- k_{0}}{Q_{L + 1}{({z - z_{L}})}}}{f_{{L + 1},l^{\prime}}^{+}\left( z^{\prime} \right)}} +} \right.} \\\left. {^{k_{0}{Q_{L + 1}{({z - z_{L}})}}}{g_{{L + 1},l^{\prime}}^{+}\left( z^{\prime} \right)}} \right\rbrack\end{matrix} & \; \\\left( {{{for}\mspace{14mu} l} = {{L + 1} > l^{\prime}}} \right) & \;\end{matrix} \right.} & {{Eq}.\mspace{14mu} (46)}\end{matrix}$

where W₀, Q₀, W_(L+1), Q_(L+1) are the same as defined previously, andf_(Ol′) ⁻(z′) and g_(L+l,l′) ⁺(z′) are zero matrices.

A. Application of Transfer Matrix Method from Substrate (l=L+1) toSource Layer (l=l′)

Applying matching conditions for G(z,z′) and G′(z, z^(l))≡∂G(z, z′)/∂zat z=z_(l) (l=l′, . . . , L), we obtain a recursion relation of thesimilar form to Eq. (36):

$\begin{matrix}{\begin{bmatrix}{f_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)} \\{g_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)}\end{bmatrix} = {{\begin{bmatrix}X_{l}^{- 1} & 0 \\0 & X_{l}\end{bmatrix}\begin{bmatrix}M_{l}^{+} & M_{l}^{-} \\M_{l}^{-} & M_{l}^{+}\end{bmatrix}}\begin{bmatrix}{f_{{l + 1},l^{\prime}}^{+}\left( z^{\prime} \right)} \\{g_{{l + 1},l^{\prime}}^{+}\left( z^{\prime} \right)}\end{bmatrix}}} & {{Eq}.\mspace{14mu} (47)}\end{matrix}$

where X_(l) and M_(l) ^(±) are as defined in Eqs. (37) and (38).Repeated applications of Eq. (47) will eventually yield:

$\begin{matrix}{\begin{bmatrix}{f_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)} \\{g_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)}\end{bmatrix} = {\begin{bmatrix}a_{l}^{+} \\b_{l}^{+}\end{bmatrix}{f_{{L + 1},l^{\prime}}^{+}\left( z^{\prime} \right)}}} & {{Eq}.\mspace{14mu} (48)}\end{matrix}$

Eliminating f_(L+1,l′) ⁺(z′), we obtain

g _(ll′) ⁺(z′)=b _(l) ⁺(a _(l) ⁺)⁻¹ f _(ll′) ⁺(z′)≡R _(l) f _(ll′)⁺(z′)  Eq. (49)

Substituting Eq. (49) into Eq. (47), we obtain the same recursionrelation as in Eq. (41)

R _(l) =X _(l)(M _(l) ⁻ +M _(l) ⁺ R _(l+1))(M _(l) ⁺ M _(l) ⁻ R _(l+1))⁻X _(l)  Eq. (50)

Because the initial setting value R_(L+1)=O is the same, we can useR_(l) determined in section II, without repeating the same calculation.

Since we have Eq. (49) relating g_(ll′) ⁺(z′) and f_(ll′) ⁺(z′) and Eq.(50) giving the algorithm to determine R_(l), f_(ll′) ⁺(z′) should besolved to obtain G_(ll′)(z,z′). As a first step, we will express f_(ll′)⁺(z′) in terms of f_(l′l′) ⁺(z′) . To this end, Eq. (47) is inverted tohave

$\begin{matrix}{\begin{bmatrix}{f_{{l + 1},l^{\prime}}^{+}\left( z^{\prime} \right)} \\{g_{{l + 1},l^{\prime}}^{+}\left( z^{\prime} \right)}\end{bmatrix} = {{\begin{bmatrix}N_{l}^{+} & N_{l}^{-} \\N_{l}^{-} & N_{l}^{+}\end{bmatrix}\begin{bmatrix}X_{l} & 0 \\0 & X_{l}^{- 1}\end{bmatrix}}\begin{bmatrix}{f_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)} \\{g_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)}\end{bmatrix}}} & {{Eq}.\mspace{14mu} (51)}\end{matrix}$

which, together with Eq. (49), gives a recursion relation for f_(ll′)⁺(z′):

$\begin{matrix}\begin{matrix}{{f_{{l + 1},l^{\prime}}^{+}\left( z^{\prime} \right)} = {\left( {{N_{l}^{+}X_{l}} + {N_{l}^{-}X_{l}R_{l}}} \right){f_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)}}} \\{\equiv {\Psi_{l}{f_{{ll}^{\prime}}^{+}\left( z^{\prime} \right)}}}\end{matrix} & {{Eq}.\mspace{14mu} (52)}\end{matrix}$

Applying Eq. (52) repeatedly, for l=l′+1, . . . , L, we obtain

f _(ll′) ⁺(z′)=Ψ_(l−1) . . . Ψ_(l′) f _(l′l′) ⁺(z′)  Eq. (53)

Once f_(l′l′) ⁺(z′) is obtained, all of f_(ll′) ⁺(z′) and g_(ll′) ⁺(z′)are determined. Before obtaining f_(l′l′) ⁺(z′), we return to theproblem of obtaining G_(ll′)(z,z′) for l<l′.

B. Application of Transfer Matrix Method from Incident Medium (l=O) toSource Layer (l=l′)

Matching conditions for G(z,z′) and G′(z,z′)≡∂G(z, z′)/∂z at z=z_(l)(l=O, . . . , l′−1) yield

$\begin{matrix}{\begin{bmatrix}{f_{{l + 1},l^{\prime}}^{-}\left( z^{\prime} \right)} \\{g_{{l + 1},l^{\prime}}^{-}\left( z^{\prime} \right)}\end{bmatrix} = {{\begin{bmatrix}N_{l}^{+} & N_{l}^{-} \\N_{l}^{-} & N_{l}^{+}\end{bmatrix}\begin{bmatrix}X_{l} & 0 \\0 & X_{l}^{- 1}\end{bmatrix}}\begin{bmatrix}{f_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)} \\{g_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)}\end{bmatrix}}} & {{Eq}.\mspace{14mu} (54)}\end{matrix}$

where X_(l) and N_(l) ^(±)are as defined previously. Using Eq. (54)repeatedly, we obtain

$\begin{matrix}{\begin{bmatrix}{f_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)} \\{g_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)}\end{bmatrix} = {\begin{bmatrix}a_{l}^{-} \\b_{l}^{-}\end{bmatrix}{g_{0l^{\prime}}^{-}\left( z^{\prime} \right)}}} & {{Eq}.\mspace{14mu} (55)}\end{matrix}$

By elimination of g_(Ol′) ⁻(z′) in this equation, f_(ll′) ⁻(z′) can beexpressed in terms of g_(ll′) ⁻(z′):

f _(ll′) ⁻(z′)=a _(l) ⁻(b _(l) ⁻)⁻¹ g _(ll′) ⁻(z′)≡ r _(l) g _(ll′)⁻(z′)  Eq. (56)

Eliminating f_(l+1,l′) ⁻(z′) and f_(ll′) ⁻(z′) in Eq. (54) by using Eq.(56), we obtain a recursion relation:

r _(l+1)=(N _(l) ⁻ +N _(l) ⁺ X _(l) r _(l) X _(l))(N _(l) ⁺ +N _(l) ⁻ X_(l) ⁻¹ r _(l) X _(l))⁻¹  Eq. (57)

By repeated applications of this recursion relation starting from l=Owith the initial setting r ₀=O, we obtain all values of r _(l).

Since f_(ll′) ⁻(z′) is expressed in terms of g_(ll′) ⁻(z′), theremaining problem is to determine g_(ll′) ³¹ (z′). To this end, g_(ll′)⁻(z′) is expressed in terms of g_(l′l′) ⁻(z′). By application of theboundary conditions at z=z_(l) (l<l′), a matching equation of the formgiven in Eq. (47) is obtained

$\begin{matrix}{\begin{bmatrix}{f_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)} \\{g_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)}\end{bmatrix} = {{\begin{bmatrix}X_{l}^{- 1} & 0 \\0 & X_{l}\end{bmatrix}\begin{bmatrix}M_{l}^{+} & M_{l}^{-} \\M_{l}^{-} & M_{l}^{+}\end{bmatrix}}\begin{bmatrix}{f_{{l + 1},l^{\prime}}^{-}\left( z^{\prime} \right)} \\{g_{{l + 1},l^{\prime}}^{-}\left( z^{\prime} \right)}\end{bmatrix}}} & {{Eq}.\mspace{14mu} (58)}\end{matrix}$

From Eq. (58) and Eq. (56), we obtain a recursion relation for g_(ll′)⁻(z′):

$\begin{matrix}\begin{matrix}{{g_{{ll}^{\prime}}^{-}\left( z^{\prime} \right)} = {{X_{l}\left( {{M_{l}^{-}{\overset{\_}{r}}_{l + 1}} + M_{l}^{+}} \right)}{g_{{l + 1},l^{\prime}}^{-}\left( z^{\prime} \right)}}} \\{\equiv {\Phi_{l + 1}{g_{{l + 1},l^{\prime}}^{-}\left( z^{\prime} \right)}}}\end{matrix} & {{Eq}.\mspace{14mu} (59)}\end{matrix}$

which enables us to express g_(ll′) ⁻(z′) (l=1, . . . , l′−1) in termsof g_(l′l′) ⁻(z′):

g _(ll′) ⁻(z′)=Φ_(l+1) . . . Ψ_(l′) g _(l′l′) ⁻(z′)  Eq. (60)

When g_(l′l′) ⁻(z′) is obtained, all values of g_(ll′) ⁻(z′) and f_(ll′)⁻(z′) (l=1, . . . , l′−1) are determined.

C. Application of Boundary Conditions of G and ∂G/∂z at z=z′

Now let us fix f_(l′l′) ⁺(z′) and g_(l′l′) ⁻. Since G_(l′l′)(z,z′) iscontinuous and G_(l′l′)′(z, z′)/∂z is the discontinuous at z=z′ by adelta function on the right-hand side in Eq. (44), we obtain

$\begin{matrix}\begin{matrix}{\begin{bmatrix}{f_{l^{\prime}l^{\prime}}^{+}\left( z^{\prime} \right)} \\{g_{l^{\prime}l^{\prime}}^{-}\left( z^{\prime} \right)}\end{bmatrix} = {\frac{k_{0}}{2}\begin{bmatrix}\left\lbrack {I - {{\overset{\_}{r}}_{l^{\prime}}R_{l^{\prime}}}} \right\rbrack^{- 1} & {- {{\overset{\_}{r}}_{l^{'}}\left\lbrack {I - {R_{l^{\prime}}{\overset{\_}{r}}_{l^{\prime}}}} \right\rbrack}^{- 1}} \\{R_{l^{\prime}}\left\lbrack {I - {{\overset{\_}{r}}_{l^{\prime}}R_{l^{\prime}}}} \right\rbrack}^{- 1} & {- \left\lbrack {I - {R_{l^{\prime}}{\overset{\_}{r}}_{l^{\prime}}}} \right\rbrack^{- 1}}\end{bmatrix}}} \\{\begin{bmatrix}^{k_{0}{Q_{l^{\prime}}{({z^{\prime} - z_{l^{\prime} - 1}})}}} \\{- ^{{- k_{0}}{Q_{l^{\prime}}{({z^{\prime} - z_{l^{\prime} - 1}})}}}}\end{bmatrix}Q_{l^{\prime}}^{- 1}W_{l^{\prime}}^{T}}\end{matrix} & {{Eq}.\mspace{14mu} (61)}\end{matrix}$

With u_(l′)=(I−R_(l′) r _(l′))⁻¹, Eq. (61) can be written as

$\begin{matrix}\begin{matrix}{{f_{l^{\prime}l^{\prime}}^{+}\left( z^{\prime} \right)} = {\frac{k_{0}}{2}\left\lbrack {^{k_{0}{Q_{l^{\prime}}{({z^{\prime} - z_{l^{\prime} - 1}})}}} + {{\overset{\_}{r}}_{l^{\prime}}{u_{l^{\prime}}\left\lbrack {{R_{l^{\prime}}^{k_{0}{Q_{l^{\prime}}{({z^{\prime} - z_{l^{\prime} - 1}})}}}} +} \right.}}} \right.}} \\{\left. \left. ^{{- k_{0}}{Q_{l^{\prime}}{({z^{\prime} - z_{l^{\prime} - 1}})}}} \right\rbrack \right\rbrack Q_{l^{\prime}}^{- 1}W_{l^{\prime}}^{T}}\end{matrix} & {{Eq}.\mspace{14mu} (62)} \\\begin{matrix}{{g_{l^{\prime}l^{\prime}}^{-}\left( z^{\prime} \right)} = {\frac{k_{0}}{2}{u_{l^{\prime}}\left\lbrack {{R_{l^{\prime}}^{k_{0}{Q_{l^{\prime}}{({z^{\prime} - z_{l^{\prime} - 1}})}}}} +} \right.}}} \\{\left. ^{{- k_{0}}{Q_{l^{\prime}}{({z^{\prime} - z_{l^{\prime} - 1}})}}} \right\rbrack Q_{l^{\prime}}^{- 1}W_{l^{\prime}}^{T}}\end{matrix} & {{Eq}.\mspace{14mu} (63)}\end{matrix}$

From Eqs. (62), (63), (53) and (60), the Green's function of Eq. (45) iscompletely determined.

IV. Numerical Implementation: TE Mode

A. Discretization of Lippmann-Schwinger Equation

Two quantities, S_(y) ^([0])(z) and G(z,z′), which are defined by thephysical properties and geometrical formation of the zeroth orderperiodic structure and the information of the incident wave, and theperturbation potential Y(z) are input parameters. An unknown quantityS_(y)(z) to be determined by the physical properties and geometricalformation of the perturbed periodic structure and the information of thesame incident wave is to be obtained by solving the Lippmann-Schwingerequation iteratively.

It is assumed that the perturbation function Y(z′) has a constant valueY_(j) inside layer j. By setting z→z_(l−1) for z locating inside thelayer l (l=1, . . . , L) and z→z_(L) for z locating inside the layer Land by using Eqs. (53) and (60), we transform Eq. (9) as follows:

$\begin{matrix}\begin{matrix}{{{S_{y}\left( z_{l - 1} \right)} - {S_{y}^{\lbrack 0\rbrack}\left( z_{l - 1} \right)}} =} \\{W_{l}\left\lbrack {{\left( {1 + R_{l}} \right){\sum\limits_{j = 1}^{l - 1}{\Psi_{l - 1}\mspace{11mu} \ldots \mspace{11mu} \Psi_{j}{\int_{z_{j - 1}}^{z_{j}}{{z^{\prime}}{f_{jj}^{+}\left( z^{\prime} \right)}Y_{j}S_{y}\left( z^{\prime} \right)}}}}} +} \right.} \\{{\left( {{\overset{\_}{r}}_{l} + 1} \right){\int_{z_{l - 1}}^{z_{l}}{{z^{\prime}}{g_{ll}^{-}\left( z^{\prime} \right)}Y_{l}S_{y}\left( z^{\prime} \right)}}} +} \\\left. {\left( {{\overset{\_}{r}}_{l} + 1} \right){\sum\limits_{j = {l + 1}}^{L}{\Phi_{l + 1}\mspace{11mu} \ldots \mspace{11mu} \Phi_{j}{\int_{z_{j - 1}}^{z_{j}}{{z^{\prime}}{g_{jj}^{-}\left( z^{\prime} \right)}Y_{j}{S_{y}\left( z^{\prime} \right)}}}}}} \right\rbrack\end{matrix} & {{Eq}.\mspace{14mu} (64)} \\\begin{matrix}{{{S_{y}\left( z_{L} \right)} - {S_{y}^{\lbrack 0\rbrack}\left( z_{L} \right)}} = {W_{L}\left( {X_{L} + {X_{L}^{- 1}R_{L}}} \right)}} \\{\quad\left\lbrack {{\sum\limits_{j = 1}^{L - 1}{\Psi_{L - 1}\mspace{11mu} \ldots \mspace{11mu} \Psi_{j}{\int_{z_{j - 1}}^{z_{j}}{{z^{\prime}}{f_{jj}^{+}\left( z^{\prime} \right)}Y_{j}S_{y}\left( z^{\prime} \right)}}}} +} \right.} \\\left. {\int_{z_{L - 1}}^{z_{L}}{{z^{\prime}}{f_{LL}^{+}\left( z^{\prime} \right)}Y_{L}{S_{ll}\left( z^{\prime} \right)}}} \right\rbrack\end{matrix} & {{Eq}.\mspace{14mu} (65)}\end{matrix}$

It is noted that, in Eq. (64), if l=1,

Σ_(j − 1)^(l − 1)(  …  )

generates no terms and if l=L,

Σ_(j = l + 1)^(L)(  …  )

generates no terms.

To enhance accuracy, we carry out analytic integrations within eachsegment by using the interpolation formula

$\begin{matrix}{{S_{y}\left( z^{\prime} \right)} = {{S_{y}\left( z_{j - 1} \right)} + {{S_{y}^{\prime}\left( z_{j} \right)}\left( {z^{\prime} - z_{j - 1}} \right)\left( {{{z_{j - 1} < z^{\prime} < z_{j}};{j = 1}},\ldots \mspace{14mu},L} \right)\mspace{14mu} {where}}}} & {{Eq}.\mspace{14mu} (66)} \\{{{{S_{y}^{\prime}\left( z_{j} \right)} \equiv {\left\lbrack {{S_{y}\left( z_{j} \right)} - {S_{y}\left( z_{j - 1} \right)}} \right\rbrack/d_{j}}},{d_{j} \equiv {z_{j} - z_{j - 1}}}}{{{{Eq}.\mspace{14mu} (64)}\mspace{14mu} {becomes}},{{{for}\mspace{14mu} l} = 2},\ldots \mspace{14mu},{L - 1},}} & {{Eq}.\mspace{14mu} (67)} \\{{{{S_{y}\left( z_{l - 1} \right)} - {S_{y}^{\lbrack 0\rbrack}\left( z_{l - 1} \right)}} = {{\sum\limits_{j = 0}^{l - 2}{\Theta_{{l - 1},j}^{( + )}Y_{j + 1}{S_{y}\left( z_{j} \right)}}} + {\sum\limits_{j = {l - 1}}^{L - 1}{\Theta_{{l - 1},j}^{( - )}Y_{j + 1}{S_{y}\left( z_{j} \right)}}} + {\sum\limits_{j = 1}^{l - 1}{{\overset{\_}{\Theta}}_{{l - 1},j}^{( + )}Y_{j}{S_{y}\left( z_{j} \right)}}} + {\sum\limits_{j = 1}^{L}{{\overset{\_}{\Theta}}_{{l - 1},j}^{( - )}Y_{j}{S_{y}\left( z_{j} \right)}}}}},\mspace{14mu} {{{for}\mspace{14mu} l} = 1},} & {{Eq}.\mspace{14mu} (68)} \\{{{{{S_{y}\left( z_{0} \right)} - {S_{y}^{\lbrack 0\rbrack}\left( z_{0} \right)}} = {{\sum\limits_{j = 0}^{L - 1}{\Theta_{0_{j}}^{( - )}Y_{j + 1}{S_{y}\left( z_{j} \right)}}} + {\sum\limits_{j = 1}^{L}{{\overset{\_}{\Theta}}_{0_{j}}^{( - )}Y_{j}{S_{y}\left( z_{j} \right)}}}}},{{{for}\mspace{14mu} l} = L},}\mspace{14mu}} & {{Eq}.\mspace{14mu} (69)} \\{{{{S_{y}\left( z_{L - 1} \right)} - {S_{y}^{\lbrack 0\rbrack}\left( z_{L - 1} \right)}} = {{\sum\limits_{j = 0}^{L - 2}{\Theta_{{L - 1},j}^{( + )}Y_{j + 1}{S_{y}\left( z_{j} \right)}}} + {\Theta_{{L - 1},{L - 1}}^{( - )}Y_{L}{S_{y}\left( z_{L - 1} \right)}} + {\sum\limits_{j = 1}^{L - 1}{{\overset{\_}{\Theta}}_{{L - 1},j}^{( + )}Y_{j}{S_{y}\left( z_{j} \right)}}} + {{\overset{\_}{\Theta}}_{{L - 1},L}^{( - )}Y_{L}{S_{y}\left( z_{L} \right)}}}}{{{Eq}.\mspace{14mu} (65)}\mspace{14mu} {becomes}}} & {{Eq}.\mspace{14mu} (70)} \\{{{S_{y}\left( z_{L} \right)} - {S_{y}^{\lbrack 0\rbrack}\left( z_{L} \right)}} = {{\sum\limits_{j = 0}^{L - 1}{\Theta_{L,j}^{( + )}Y_{j + 1}{S_{y}\left( z_{j} \right)}}} + {\sum\limits_{j = 1}^{L}{{\overset{\_}{\Theta}}_{L,j}^{( + )}Y_{j}{S_{y}\left( z_{j} \right)}}}}} & {{Eq}.\mspace{14mu} (71)}\end{matrix}$

Eqs. (68), (69), (70) and (71) can be put together in one expandedmatrix form of a linear system of equations:

=[

−

+

  Eq. (72)

where

and

are column vectors of dimension L+1 with the layer componentsS_(y)(z_(l)) and S_(y) ^([0])(z_(l)) respectively. Each layer componentS_(y)(z_(l)) and S_(y) ^([0])(z_(l)) are still column vectors ofdimension (2N+1) with the coupled-wave basis components.

are square matrices with (L+1)² layer components, each layer componentbeing (2N+1)² square matrices in coupled-wave:

$\begin{matrix}{= \left\{ \begin{matrix}\Theta_{lm}^{( + )} & {{{{for}\mspace{14mu} l} > m},{m \neq L}} \\\Theta_{lm}^{( - )} & {{{{for}\mspace{14mu} l} \leq m},{m \neq L}} \\0 & {{{for}\mspace{14mu} m} = L}\end{matrix} \right.} & {{Eq}.\mspace{14mu} (73)} \\{= \left\{ \begin{matrix}0 & {{{for}\mspace{14mu} l} - 0} \\{\overset{\_}{\Theta}}_{lm}^{( + )} & {{{{for}\mspace{14mu} l} \neq 0},{l \geq m}} \\\Theta_{lm}^{( - )} & {{{{for}\mspace{14mu} l} \neq 0},{l < m}}\end{matrix} \right.} & {{Eq}.\mspace{14mu} (74)} \\{= \left\{ \begin{matrix}{Y_{l + 1}\delta_{lm}} & {{{for}\mspace{14mu} l} \neq L} \\0 & {{{for}\mspace{14mu} l} = L}\end{matrix} \right.} & {{Eq}.\mspace{14mu} (75)} \\{= \left\{ \begin{matrix}0 & {{{for}\mspace{14mu} l} = 0} \\{Y_{l}\delta_{lm}} & {{{for}\mspace{14mu} l} \neq 0}\end{matrix} \right.} & {{Eq}.\mspace{14mu} (76)}\end{matrix}$

B. Special Cases

As special cases,

-   -   (i) If each layer is a homogeneous medium (case 1), we obtain

W₁= . . . =W_(L)=I  Eq. (77)

whereby R_(l), r _(l), u_(l), Ψ_(l), Φ_(l) become diagonal matrixes. Asa result, since matrix elements

and

in Eqs. (73) and (74) all become diagonal matrixes, we obtain moresimplified formulas.

-   -   (ii) If a zeroth order dielectric function of each layer is the        same (case 4), that is, when

W₁= . . . =W_(L)=W  Eq. (78)

Q₁= . . . =Q_(L)=Q  Eq. (79)

whereby simplified matrix elements

and

can be obtained.

(iii) Furthermore, if each layer is homogeneous (case 2),

W₁= . . . =W_(L)=I  Eq. (80)

Q₁= . . . =Q_(L)=qI  Eq. (81)

whereby more simply expressed matrix elements can be obtained.

-   -   (iv) If the perturbation region is taken as the incident medium        (case 3), this corresponds to case 2 with u=O. Thus, the matrix        elements become much simpler.

C. Calculation of S_(y)(z)

To enhance the speed of calculation of multiplying a matrix and a vectorin the space of the coupled-wave basis, the perturbation matrix can beblock-diagonalized from a similarity transformation with

$\begin{matrix}{T = {\frac{1}{2}\begin{bmatrix}{- I_{N}} & 0_{c} & J_{N} \\0_{r} & I_{1} & 0_{r} \\J_{N} & 0_{c} & I_{N}\end{bmatrix}}} & {{Eq}.\mspace{14mu} (82)}\end{matrix}$

where I_(K) is a K×K unit matrix, J_(N) is an N×N square matrix withelements J_(ij)=δ_(i,N+1−j), and 0_(c) and 0_(r) are, respectively, zerocolumn and zero row vectors of dimension N.

To solve the matrix equation of Eq. (72) iteratively, we start with aninitial settings

=

. We substitute this into Eq. (72),

=

−[

+

  Eq. (83)

If both sides of this resulting equation are not equal within a givenerror range, we substitute a new

defined as the right-hand side of Eq. (83),

=

−[

+

  Eq. (84)

into Eq. (72) to obtain

−[

+

=

−[

+

](

−

+

)  Eq. (85)

We compare both sides of this equation. If the comparison does not givea satisfactory result, we define a new

as the right-hand side of Eq. (85):

=

−[

+

](

−

+

)  Eq. (86)

With this

we test the validity of Eq. (72). This process of comparison will becontinued until a satisfactory result is obtained.

D. Calculation of g₀ and f_(L+1)

The final goal we wish to achieve is to calculate reflectivity ρ_(n) andtransmittance τ_(n) related to the order n of the coupled-wave basis.The reflectivity column vector g₀=(ρ_(−N), . . . , ρ⁻¹, ρ₀, ρ₁, . . . ,ρ_(N))^(T) can be extracted as follows:

$\begin{matrix}{g_{0} = {g_{0}^{\lbrack 0\rbrack} + {\sum\limits_{l = 0}^{L}{\sum\limits_{m = 0}^{L}\left\lbrack + \right.}}}} & {{Eq}.\mspace{14mu} (87)}\end{matrix}$

g₀ ^([0])=(ρ_(−N) ^([0]), . . . , ρ⁻¹ ^([0]), ρ₀ ^([0]), ρ₁ ^([0]), . .. , ρ_(N) ^([0]))^(T) is the zeroth order reflectivity column vector.The middle component of Eq. (87), ρ₀=(g₀)₀, characterizes a reflectiondue to the principal order.

Similarly, the transmittance column vector f_(L+1)=(τ_(−N), . . . , τ⁻¹,τ₀, τ₁, . . . , τ_(N))^(T) is given as follows:

$\begin{matrix}{f_{L + 1} = {f_{L + 1}^{\lbrack 0\rbrack} + {\sum\limits_{l = 0}^{L}{\sum\limits_{m = 0}^{L}\left\lbrack + \right.}}}} & {{Eq}.\mspace{14mu} (88)}\end{matrix}$

where the components of zeroth order transmittance column vector aref_(L+1) ^([0])=(τ_(−N) ^([0]), . . . , τ⁻¹ ^([0]), τ₀ ^([0]), τ₁ ^([0]),. . . , τ_(N) ^([0]))^(T), and τ₀=(f_(L+1))₀ is the transmittance of atransmitted wave by the principal order.

V. Lippmann-Schwinger Equation: TM Mode

Due to the periodicity of the periodic structure in x direction, the TMmode solutions to Maxwell's Eqs. can be expressed as

$\begin{matrix}{{H_{y}\left( {x,z} \right)} = {\sum\limits_{n = {- \infty}}^{\infty}{{U_{yn}(z)}^{\; k_{xn}x}}}} & {{Eq}.\mspace{14mu} (89)} \\{{E_{x}\left( {x,z} \right)} = {{- {\left( \frac{\mu_{0}}{\varepsilon_{0}} \right)}^{1/2}}{\sum\limits_{n = {- \infty}}^{\infty}{{S_{xn}(z)}^{\; k_{xn}x}}}}} & {{Eq}.\mspace{14mu} (90)}\end{matrix}$

The inverse of the dielectric function ε(x,z) can be expanded in Fourierseries as follows:

$\begin{matrix}{\frac{1}{ɛ\left( {x,z} \right)}{\sum\limits_{h = {- \infty}}^{\infty}{{ɛ_{h}^{\#}(z)}^{\; 2\pi \; {{hx}/\Lambda}}}}} & {{Eq}.\mspace{14mu} (91)}\end{matrix}$

From Maxwell's equations, U_(yn)(z) and S_(xn)(z) are coupled via twolinear differential equations. Eliminating S_(xn)(z), one can get asecond-order decoupled differential equation of U_(yn)(z):

$\begin{matrix}\begin{matrix}\left\lbrack {{\frac{1}{k_{0}^{2}}\frac{\partial^{2}}{\partial z^{2}}} - {\frac{1}{k_{0}^{2}}\frac{\partial{E(z)}}{\partial_{z}}{E^{- 1}(z)}\frac{\partial{U_{y}(z)}}{\partial z}} -} \right. \\{{\left. {E(z){B(z)}} \right\rbrack {U_{y}(z)}} = {0\mspace{14mu} {where}}}\end{matrix} & {{Eq}.\mspace{14mu} (92)} \\{{B(z)} = {{K_{x}{P(z)}K_{x}} - I}} & {{Eq}.\mspace{14mu} (93)}\end{matrix}$

where P(z) is a square matrix with the (n,p) element being equal toP_(np)(z)=ε_((n−p)) ^(#)(z). A similar equation for the zeroth orderperiodic structure with the same period can be obtained by using thesuperscript [0]:

$\begin{matrix}\begin{matrix}\left\lbrack {{\frac{1}{k_{0}^{2}}\frac{\partial^{2}}{\partial z^{2}}} - {\frac{1}{k_{0}^{2}}\frac{\partial{E^{\lbrack 0\rbrack}(z)}}{\partial z}\left( E^{\lbrack 0\rbrack} \right)^{- 1}(z)\frac{\partial{U_{y}^{\lbrack 0\rbrack}(z)}}{\partial z}} -} \right. \\{{\left. {{E^{\lbrack 0\rbrack}(z)}{B^{\lbrack 0\rbrack}(z)}} \right\rbrack {U_{y}^{\lbrack 0\rbrack}(z)}} = 0}\end{matrix} & {{Eq}.\mspace{14mu} (94)}\end{matrix}$

Corresponding Green's function is given as follows:

$\begin{matrix}\begin{matrix}\left\lbrack {{\frac{1}{k_{0}^{2}}\frac{\partial^{2}}{\partial z^{2}}} - {\frac{1}{k_{0}^{2}}\frac{\partial{E^{\lbrack 0\rbrack}(z)}}{\partial z}\left( E^{\lbrack 0\rbrack} \right)^{- 1}(z)\frac{\partial{U_{y}^{\lbrack 0\rbrack}(z)}}{\partial z}} -} \right. \\{{\left. {{E^{\lbrack 0\rbrack}(z)}{B^{\lbrack 0\rbrack}(z)}} \right\rbrack {G\left( {z,z^{\prime}} \right)}} = {{- {\delta \left( {z - z^{\prime}} \right)}}I}}\end{matrix} & {{Eq}.\mspace{14mu} (95)}\end{matrix}$

Based on the standard method, combining Eqs. (92), (94) and (95), weobtain the Lippmann-Schwinger equation for the TM mode as follows:

$\begin{matrix}\begin{matrix}{{{U_{y}(z)} - {U_{y}^{(0)}(z)}} = {\int{{z^{\prime}\left\lbrack {{k_{0}^{2}{G\left( z^{\prime} \right)}K_{x}{\overset{\sim}{Y}\left( z^{\prime} \right)}K_{x}U_{y}\left( z^{\prime} \right)} +} \right.}}}} \\{{\left. {\frac{\partial{G\left( {z,z^{\prime}} \right)}}{\partial z^{\prime}}{Y\left( z^{\prime} \right)}\frac{\partial{U_{y}\left( z^{\prime} \right)}}{\partial z^{\prime}}} \right\rbrack {Y(z)}} \equiv {{E^{- 1}(z)} - {\left( E^{\lbrack 0\rbrack} \right)^{- 1}(z)}}}\end{matrix} & {{Eq}.\mspace{14mu} (96)}\end{matrix}$

where

VI. Calculation of U_(y) ^([0])(z)

In this section, we omit the superscript [0] used in the previoussection to indicate the relevant quantities in the zeroth order periodicstructure.

Solutions to the electric field and the magnetic field in region I andregion II can be expanded in the coupled-wave basis by the Floquet'sconditions:

$\begin{matrix}{{H_{y}^{(0)}\left( {x,z} \right)} = {\sum\limits_{n = {- N}}^{N}{{U_{yn}^{(0)}(z)}^{\; k_{xn}x}}}} & {{Eq}.\mspace{14mu} (97)} \\{{E_{x}^{(0)}\left( {x,z} \right)} = {{- {\left( \frac{\mu_{0}}{\varepsilon_{0}} \right)}^{1/2}}{\sum\limits_{n = {- N}}^{N}{{S_{xn}^{(0)}(z)}^{\; k_{zn}x}}}}} & {{Eq}.\mspace{14mu} (98)} \\{{H_{y}^{({L + 1})}\left( {x,z} \right)} = {\sum\limits_{n = {- N}}^{N}{{U_{yn}^{({L + 1})}(z)}^{\; k_{xn}x}}}} & {{Eq}.\mspace{14mu} (99)} \\{{E_{x}^{({L + 1})}\left( {x,z} \right)} = {{- {\left( \frac{\mu_{0}}{\varepsilon_{0}} \right)}^{1/2}}{\sum\limits_{n = {- N}}^{N}{{S_{xn}^{({L + 1})}(z)}^{\; k_{xn}x}\mspace{14mu} {where}}}}} & {{Eq}.\mspace{14mu} (100)} \\{{U_{yn}^{(0)}(z)} = {{\delta_{n0}^{\; k_{1,{zn}}z}} + {{\rho \;}_{n}^{{- }\; k_{1,{zn}}z}}}} & {{Eq}.\mspace{14mu} (101)} \\{{S_{xn}^{(0)}(z)} = {{{\delta_{n0}\left( \frac{\; k_{1,{zn}}}{n_{1}^{2}k_{0}} \right)}^{\; k_{1,{zn}}z}} - {{\rho_{n}\left( \frac{{k}_{1,{zn}}}{n_{1}^{2}k_{0}} \right)}^{{- }\; k_{1,{zn}}z}}}} & {{Eq}.\mspace{14mu} (102)} \\{{U_{yn}^{({L + 1})}(z)} = {T_{n}^{\; {k_{11,{zn}}{({z - z_{L}})}}}}} & {{Eq}.\mspace{14mu} (103)} \\{{S_{xn}^{({L + 1})}(z)} = {{T_{n}\left( \frac{\; k_{11,{zn}}}{n_{11}^{2}k_{0}} \right)}^{\; {k_{11,{zn}}{({z - z_{L}})}}}}} & {{Eq}.\mspace{14mu} (104)}\end{matrix}$

The first term on the right-hand side of Eq. (101) represents theincident wave. Likewise, solutions of the electric field and themagnetic field in the periodic structure region (l=1, . . . , L) can beexpanded in terms of coupled-wave bases as follows:

$\begin{matrix}{{H_{y}^{(l)}\left( {x,z} \right)} = {\sum\limits_{n = {- N}}^{N}{U_{yn}^{(l)}\left( {z\; } \right)}^{\; k_{xn}x}}} & {{Eq}.\mspace{14mu} (105)} \\{{E_{x}^{(l)}\left( {x,z} \right)} = {{- {\left( \frac{\mu_{0}}{\varepsilon_{0}} \right)}^{1/2}}{\sum\limits_{n = {- N}}^{N}{{S_{xn}^{(l)}(z)}^{\; k_{xn}x}}}}} & {{Eq}.\mspace{14mu} (106)}\end{matrix}$

Applying Eq. (92) to the layer l, we obtain an equation for U_(yn)^((l)) (n=−N, . . . , N) in Eq. (105):

(1/k ₀ ²)∂² U _(y) ^((l)) /∂z ² =E _(l) B _(l) U _(y) ^((l))  Eq. (107)

where E_(l) is as defined in section II, and B₁, being independent ofz_(l) is defined as follows:

B _(l) =K _(x) E _(l) ⁻¹ K _(x) −I  Eq. (108)

Solutions to Eq. (107) are given as

U _(y) ^((l))(z)=W _(l) [e ^(−k) ⁰ ^(Q) ^(l) ^((z−z) ^(l−1) ⁾ f _(l) +e^(k) ⁰ ^(Q) ^(l) ^((z−z) ^(l−1) ⁾ g _(l)]  Eq. (109)

S _(x) ^((l))(z)=W _(l) [e ^(−k) ⁰ ^(Q) ^(l) ^((z−z) ^(l−1) ⁾ f _(l) +e^(k) ⁰ ^(Q) ^(l) ^((z−z) ^(l−1) ⁾ g _(l)]  Eq. (110)

where W_(l) is a square matrix constructed from (2N+1) eigenvectors ofE_(l)B_(l), Q_(l) is a diagonal matrix with the diagonal elements beingthe positive square roots of (2N+1) eigenvalues of E_(l)B_(l),V_(l)=E_(l) ⁻¹W_(l)Q_(l), and f _(l) and g_(l) are column vectors with(2N+1) components to be determined by application of the boundaryconditions.

The components of the coupled-wave basis in Eqs. (101), (102), (103) and(104), respectively, can be expressed in the matrix form as follows:

U _(y) ⁽⁰⁾(z)=W ₀ [e ^(−k) ⁰ ^(Q) ⁰ ^(z) f ₀ +e ^(k) ⁰ ^(Q) ⁰ ^(z) g₀]  Eq. (111)

S _(x) ⁽⁰⁾(z)=V ₀ [e ^(−k) ⁰ ^(Q) ⁰ ^(z) f ₀ +e ^(k) ⁰ ^(Q) ⁰ ^(z) g₀]  Eq. (112)

U _(y) ^((L+1))(z)=W _(L+1) [e ^(−k) ⁰ ^(Q) ^(L+1) ^((z−z) ^(L) ⁾ f_(L+1) +e ^(k) ⁰ ^(Q) ^(L+1) ^((z−z) ^(L) ⁾ g _(L+1)]  Eq. (113)

S _(x) ^((L+1))(z)=V _(L+1) [e ^(−k) ⁰ ^(Q) ^(L+1) ^((z−z) ^(L) ⁾ f_(L+1) +e ^(k) ⁰ ^(Q) ^(L+1) ^((z−z) ^(L) ⁾ g _(l)]  Eq. (114)

where Q₀, W₀, f₀, g₀, W_(L+1), Q_(L+1), f_(L+1) and g_(L+1) are asdefined in section II, V₀=(1/n_(I) ²)W₀Q₀ and V_(l+1)=(1/n_(II)²)W_(L+1)Q_(L+1).

Applications of boundary conditions for the electric field and magneticfield at z=z_(l) (l=0, . . . , L) yield a recursion relation given inEq. (36). Rest part of the calculation can be carried out in the exactlysame manner as in the TE mode. We end this section with a comment on theinversion matrix of W_(l) appearing in M_(l) ^(±).

Although E_(l) and B_(l) are symmetric matrixes, the product E_(l)B_(l)is not a symmetric one in general. Thus, the eigenvectors w_(m) ^((l))(m=−N, . . . , N) of E_(l)B_(l) are not orthogonal, i.e., [w_(m)^((l))]^(T)w_(n) ^((l))≠0 for m≠n. Thus, in this case, W_(l) ^(T)W_(l)cannot be expressed as a unit matrix I by the proper normalization.Instead, since W_(l) ^(T)E_(l) ⁻¹W_(l) can be a diagonal matrix, it canbe normalized as W_(l) ^(T)E_(l) ⁻¹W_(l)=I. In TM mode, W_(l) ^(T)E_(l)⁻¹ instead of W_(l) ⁻¹ may be used.

VII. Calculation of G(z,z′): TM Mode

If z locates inside the layer l, Eq. (95) becomes

$\begin{matrix}{{\left\lbrack {{\frac{1}{k_{0}^{2}}\frac{\partial^{2}}{\partial z^{2}}} - {E_{l}B_{l}}} \right\rbrack {G\left( {z - z^{\prime}} \right)}} = {{- {\delta \left( {z - z^{\prime}} \right)}}I}} & {{Eq}.\mspace{14mu} (115)}\end{matrix}$

In this section, we omit, as in sections II and VI, the superscript [0]used to indicate the relevant quantities in the zeroth order periodicstructure. For a given layer l, E_(l) and B_(l) are not functions of zany more. Except for W_(l) and Q_(l) being the matrices constructed fromeigenvectors and eigenvalues of E_(l)B_(l), the solutions to Eq. (115)are of the same form as Eq. (45) in the TE mode, with the sameexpressions for f_(ll′) ^(±)(z′) and g_(ll′) ^(±)(z′). In this section,W_(l) and Q_(l) are exactly the same as those obtained in section VI.The solutions in the incident medium denoted as region I and thesubstrate denoted as region II are of the same form as defined in Eq.(46).

VIII. Numerical Implementation: TM Mode

A. Calculation of U_(y)(z):

To enhance accuracy of the numerical calculation, when an interpolationformula for U_(y)(z′) similar to Eq. (66) is used and analyticintegrations are carried out, the discretization of theLippmann-Schwinger equation [Eq. (97)] can be done as follows:

$\begin{matrix}\begin{matrix}{{{U_{y}\left( z_{l - 1} \right)} - {U_{y}^{\lbrack 0\rbrack}\left( z_{l - 1} \right)}} =} \\{{\sum\limits_{j = 0}^{l - 2}{\left\{ {{\Theta_{{l - 1},j}^{( + )}K_{x}{\overset{\sim}{Y}}_{j + 1}K_{x}} + {\Delta_{{l - 1},j}^{( + )}Y_{j + 1}}} \right\} U_{y}\left( z_{j} \right)}} +} \\{{\sum\limits_{j = {l - 1}}^{L - 1}{\left\{ {{\Theta_{{l - 1},j}^{( - )}K_{x}{\overset{\sim}{Y}}_{j + 1}K_{x}} + {\Delta_{{l - 1},j}^{( - )}Y_{j + 1}}} \right\} U_{y}\left( z_{j} \right)}} +} \\{{\sum\limits_{j = 1}^{l - 1}{\left\{ {{{\overset{\_}{\Theta}}_{{l - 1},j}^{( + )}K_{x}{\overset{\sim}{Y}}_{j}K_{x}} + {{\overset{\_}{\Delta}}_{{l - 1},j}^{( + )}Y_{j}}} \right\} U_{y}\left( z_{j} \right)}} +} \\\begin{matrix}{\sum\limits_{j = 1}^{L}\left\{ {{{\overset{\_}{\Theta}}_{{l - 1},j}^{( - )}K_{x}{\overset{\sim}{Y}}_{j}K_{x}} + {{\overset{\_}{\Delta}}_{{l - 1},j}^{( - )}Y_{j}}} \right\}} \\{{{{U_{y}\left( z_{j} \right)}\mspace{14mu} {for}\mspace{14mu} l} = 2},\ldots \mspace{14mu},{L - 1}}\end{matrix}\end{matrix} & {{Eq}.\mspace{14mu} (116)} \\\begin{matrix}\begin{matrix}{{{U_{y}\left( z_{0} \right)} - {U_{y}^{\lbrack 0\rbrack}\left( z_{0} \right)}} =} \\{{\sum\limits_{j = 0}^{L - 1}{\left\{ {{\Theta_{0j}^{( - )}K_{x}{\overset{\sim}{Y}}_{j + 1}K_{x}} + {\Delta_{0j}^{( - )}Y_{j + 1}}} \right\} U_{y}\left( z_{j} \right)}} +}\end{matrix} \\{\sum\limits_{j = 1}^{L}{\left\{ {{{\overset{\_}{\Theta}}_{0j}^{( - )}K_{x}{\overset{\sim}{Y}}_{j}K_{x}} + {{\overset{\_}{\Delta}}_{0j}^{( - )}Y_{j}}} \right\} {U_{y}\left( z_{j} \right)}}}\end{matrix} & {{Eq}.\mspace{14mu} (117)} \\\begin{matrix}{{{U_{y}\left( z_{L - 1} \right)} - {U_{y}^{\lbrack 0\rbrack}\left( z_{L - 1} \right)}} =} \\\begin{matrix}\begin{matrix}\begin{matrix}{{\sum\limits_{j = 0}^{L - 2}{\left\{ {{\Theta_{{L - 1},j}^{( + )}K_{x}{\overset{\sim}{Y}}_{j + 1}K_{x}} + {\Delta_{{L - 1},j}^{( + )}Y_{j + 1}}} \right\} U_{y}\left( z_{j} \right)}} +} \\{{\left\{ {{\Theta_{{L - 1},{L - 1}}^{( - )}K_{x}{\overset{\sim}{Y}}_{L}K_{x}} + {\Delta_{{L - 1},{L - 1}}^{( - )}Y_{L}}} \right\} {U_{y}\left( z_{L - 1} \right)}} +}\end{matrix} \\{{\sum\limits_{j = 1}^{L - 1}{\left\{ {{{\overset{\_}{\Theta}}_{{L - 1},j}^{( + )}K_{x}{\overset{\sim}{Y}}_{j}K_{x}} + {{\overset{\_}{\Delta}}_{{L - 1},j}^{( + )}Y_{j}}} \right\} U_{y}\left( z_{j} \right)}} +}\end{matrix} \\{\left\{ {{{\overset{\_}{\Theta}}_{{L - 1},L}^{( - )}K_{x}{\overset{\sim}{Y}}_{L}K_{x}} + {{\overset{\_}{\Delta}}_{{L - 1},L}^{( - )}Y_{L}}} \right\} {U_{y}\left( z_{L} \right)}}\end{matrix}\end{matrix} & {{Eq}.\mspace{14mu} (118)} \\\begin{matrix}\begin{matrix}{{{U_{y}\left( z_{L} \right)} - {U_{y}^{\lbrack 0\rbrack}\left( z_{L} \right)}} =} \\{{\sum\limits_{j = 0}^{L - 1}{\left\{ {{\Theta_{L,j}^{( + )}K_{x}{\overset{\sim}{Y}}_{j + 1}K_{x}} + {\Delta_{L,j}^{( + )}Y_{j + 1}}} \right\} U_{y}\left( z_{j} \right)}} +}\end{matrix} \\{\sum\limits_{j = 1}^{L}{\left\{ {{{\overset{\_}{\Theta}}_{L,j}^{( + )}K_{x}{\overset{\sim}{Y}}_{j}K_{x}} + {{\overset{\_}{\Delta}}_{L,j}^{( + )}Y_{j}}} \right\} {U_{y}\left( z_{j} \right)}}}\end{matrix} & {{Eq}.\mspace{14mu} (119)}\end{matrix}$

Eqs. (116), (117), (118) and (119) can be put together in an expandedmatrix form of a system of linear equations:

−

=[

+

+

+

  Eq. (120)

where U, U^([0]) are column vectors of dimension L+1 with layercomponents U_(y)(z_(l)) and U_(y) ^([0])(z_(l)), respectively, beingstill column vectors of dimension (2N+1) with coupled-wave basiscomponents.

are as defined in Eqs. (73), (74), (75) and (76).

are square matrices with (L+1)² layer components, each layer componentbeing square matrices with (2N+1)² components in coupled-wave basis:

$\begin{matrix}{= \left\{ \begin{matrix}\Delta_{lm}^{( + )} & {{{{for}\mspace{14mu} l} > m},{m \neq L}} \\\Delta_{lm}^{( - )} & {{{{for}\mspace{14mu} l} \leq m},{m \neq L}} \\0 & {{{for}\mspace{14mu} m} = L}\end{matrix} \right.} & {{Eq}.\mspace{14mu} (121)} \\{= \left\{ \begin{matrix}0 & {{{for}\mspace{14mu} l} = 0} \\{\overset{\_}{\Delta}}_{lm}^{( + )} & {{{{for}\mspace{14mu} l} \neq 0},{l \geq m}} \\\Delta_{lm}^{( - )} & {{{{for}\mspace{14mu} l} \neq 0},{l < m}}\end{matrix} \right.} & {{Eq}.\mspace{14mu} (122)} \\{= \left\{ \begin{matrix}{K_{x}{\overset{\sim}{Y}}_{l + 1}K_{x}\delta_{lm}} & {{{for}\mspace{14mu} l} \neq L} \\0 & {{{for}\mspace{14mu} l} = L}\end{matrix} \right.} & {{Eq}.\mspace{14mu} (123)} \\{= \left\{ \begin{matrix}0 & {{{for}\mspace{14mu} l} = 0} \\{K_{x}{\overset{\sim}{Y}}_{l}K_{x}\delta_{lm}} & {{{for}\mspace{14mu} l} = {\neq 0}}\end{matrix} \right.} & {{Eq}.\mspace{14mu} (124)}\end{matrix}$

Similarly to the TE mode, we can obtain corresponding formulas for case1, case 2, case 3 and case 4.

The matrix equation of Eq. (120) can be solved by using the iterationmethod with an initial setting

=

as in the case of the TE mode.

B. Calculation of g₀ and f_(L+1)

The reflectivity column vector g₀=(ρ_(−N), . . . , ρ⁻¹, ρ₀, ρ₁, . . . ,ρ_(N))^(T) can be extracted as follows:

$\begin{matrix}\begin{matrix}{g_{0} = {g_{0}^{\lbrack 0\rbrack} +}} \\{\sum\limits_{l = 0}^{L}{\sum\limits_{m = 0}^{L}\left\lbrack {+ + +} \right.}}\end{matrix} & {{Eq}.\mspace{14mu} (125)}\end{matrix}$

where g₀ ^([0])(ρ_(−N) ^([0]), . . . , ρ⁻¹ ^([0]), ρ₀ ^([0]), ρ₁ ^([0]),. . . , ρ_(N) ^([0]))^(T) is the zeroth order reflectivity columnvector.

The reflectivity of a principal order reflection can be obtained byextracting ρ₀=(g₀)₀,

Similarly, the transmittance column vector f_(L+1)=(τ_(−N), . . . , τ⁻¹,τ₀, τ_(1,), . . . , τ_(N))^(T) is given as:

$\begin{matrix}\begin{matrix}{f_{L + 1} = {f_{L + 1}^{\lbrack 0\rbrack} +}} \\{\sum\limits_{l = 0}^{L}{\sum\limits_{m = 0}^{L}\left\lbrack {+ + +} \right.}}\end{matrix} & {{Eq}.\mspace{14mu} (126)}\end{matrix}$

where f_(L+1) ^([0])=(τ_(−N) ^([0]), . . . , τ⁻¹ ^([0]), τ₀ ^([0]), τ₁^([0]), . . . , τ_(N) ^([0]))^(T) is the zeroth order transmittancecolumn vector, and τ_(O)=(f_(L+1))_(O) is the transmittance of theprincipal order transmission.

The above-described calculated reflectivity and transmittance of each ofthe TE mode and TM mode are compared with the measured reflectivity andtransmittance. The comparison results can be applied to thenondestructive analysis of various periodic structures, for example,holographic grating structures, surface relief and multilayer gratingstructures, plane dielectric or absorptive holographic gratingstructures, random sectional dielectric and absorptive surface reliefgrating structures, two-dimensional surface relief grating structures,or anisotropic grating structures.

The invention has been described using preferred exemplary embodiments.However, it is to be understood that the scope of the invention is notlimited to the disclosed embodiments. On the contrary, the scope of theinvention is intended to include various modifications and alternativearrangements within the capabilities of persons skilled in the art usingpresently known or future technologies and equivalents. The scope of theclaims, therefore, should be accorded the broadest interpretation so asto encompass all such modifications and similar arrangements.

1. A method of determining physical properties of a multilayeredperiodic structure, comprising steps of: (a) illuminating a realperiodic structure and performing at least one of (i) measuring at leastone of reflectivity or transmittance of the real periodic structure inresponse to the illumination or (ii) measuring at least one physicalproperty related to reflectivity or transmittance of the real periodicstructure; (b) performing at least one of (i) calculating at least oneof reflectivity or transmittance of a virtual periodic structurecorresponding to the real periodic structure in response to theillumination or (ii) calculating at least one physical property relatedto at least one of reflectivity or transmittance of the virtual periodicstructure in response to the illumination, by setting the virtualperiodic structure having a repeated shape, one-dimensionally,two-dimensionally or three-dimensionally and at least a horizontallyrepeating period and a plurality of vertically stacked layers includingmiddle layers among the vertically stacked layers, wherein at leastthree substances occur in each horizontally repeated period, and bycalculating physical properties related to the reflectivity and thetransmittance of the virtual periodic structure by using a refractiveindex of each of the at least three substances forming the virtualperiodic structure; and (c) comparing the at least one physical propertyrelated to the reflectivity and the transmittance being measured in thestep (a) with the corresponding at least one physical property relatedto the at least one of reflectivity or transmittance being calculated inthe step (b).
 2. The method of claim 1, wherein the step (b) isperformed by steps including: dividing the virtual periodic structureinto a zeroth order periodic structure and a perturbed virtual periodicstructure obtained by geometrically or physically changing the zerothorder periodic structure in a perturbation region; calculating physicalproperties related to at least one of reflectivity or transmittance oflight incident on the zeroth order periodic structure; obtaining thephysical properties related to at least one of the reflectivity or thetransmittance of the perturbed virtual periodic structure by using anapproximation algorithm together with the calculated physical propertiesrelated to at least one of reflectivity or transmittance.
 3. The methodof claim 1, wherein the physical properties related to the at least oneof the reflectivity or the transmittance being measured in the step (a)are physical properties related to at least one of amplitude or phase ofa reflected wave or a transmitted wave in a TE mode electric field andphysical properties related to amplitude or phase of a reflected wave ora transmitted wave in a TM mode magnetic field, and the at least one ofthe reflectivity or the transmittance being calculated in the step (b)are the at least one of the amplitude or phase of the reflected wave orthe transmitted wave in the TE mode electric field and the amplitude orphase of the reflected wave or the transmitted wave in the TM modemagnetic field.
 4. The method of claim 1, wherein the step (b) comprisessteps of: expanding a dielectric function in the virtual periodicstructure in Fourier series; expanding the incident wave, reflected waveand transmitted wave of the beam being incident on the virtual periodicstructure in a sum of a plane wave of an electromagnetic wave; andcalculating the reflectivity and the transmittance by using boundaryconditions of expansion coefficients.
 5. The method of claim 1, whereinthe virtual periodic structure comprises a ridge region formed betweengroove regions being formed of a third substance, and the ridge regioncomprises a center part formed of a first substance and a surface layerformed of a second substance on the surface of the center part.
 6. Themethod of claim 5, wherein the surface layer includes at least one of alayer selected from the group consisting of an oxide layer, a coatinglayer, or a surface roughness layer.
 7. The method of claim 5, whereinthe third substance has at least one of a gaseous, liquid or solidphase.
 8. The method of claim 5, wherein the step (b) comprisesexpanding the dielectric function of each vertically stacked layer ofthe virtual periodic structure in the Fourier series, and the expansioncoefficient of the Fourier series in each such layer of the virtualperiodic structure is expressed by using a complex refractive index ofeach of the first substance and second substance forming the ridgeregion and the third substance forming the groove region, a ratio f_(l)of the region occupied by the first and second substances of each layerl with respect to a period Λ, a ratio 2δ_(l)/Λ of the region occupied bythe second substance only, and a parameter t_(l) indicating how far thecenter of the layer l is off from the center of layer 1 of the dividedlayers in x direction.
 9. The method of claim 5, wherein the firstsubstance is formed of at least two or more different substances beinghorizontally or vertically positioned.